Melting Transition of Vortex Lattice in Point Vortex Systems 
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Point vortices take a triangular lattice structure in a rotating system as a minimum energy state. 
We perform a numerical simulation of point vortex systems using initial conditions indicating that 
the triangular lattice is randomly perturbed. The total energy increases with the magnitude of 
the perturbation. When the energy is increased, the vortex lattice becomes irregular and a layered 
structure appears. When the energy is further increased, the layered structure disappears and a 
liquidlike state appears. We interpret the melting transition with a mean- field approximation for 
O-i layered structures. 
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I. INTRODUCTION 

Quantized vortices are important in several quantum systems such as superfluid 4 He and superconductors. Tri- 
angular lattices of vortices and magnetic fluxes were observed in rotating superfluids and type II superconductors 
in magnetic fields. Similar quantized vortex lattices have recently been found in rapidly rotating Bose-Einstein con- 
densates of ultracold diluted trapped gases. [HQ The quantized vortex states in the Bose-Einstein condensates have 
been intensively studied theoretically. [3- 5] The vortex lattice state is considered to be the lowest energy state in the 
rotating Bose-Einstein condensates with repulsive interaction. The Gross-Pitaevskii equation describes fairly well the 
Bose-Einstein condensates. The Gross-Pitaevskii equation is an equation in which the total energy and the angular 
momentum are conserved. If the initial conditions deviate from the vortex lattice state, various types of complex dy- 
namics are expected to occur. We studied the time evolution of the Gross-Pitaevskii equation from an initial condition 
analogous to a rigid rotation, and found that a vortex lattice state appeared owing to Kelvin-Helmholtz instability. [(| 
Since the obtained state is not the lowest-energy state, the vortices exhibited an irregular rotating motion around the 
regular lattice points. When the irregular motion becomes large, the lattice structure is expected to melt into a liquid 
state. 

Several phenomena similar to the melting phenomena of vortex lattices have been studied in other experiments. A 
liquid-to-crystal transition was studied in a two-dimensional sheet of electrons [7[ and a transition from turbulence 
to vortex crystals was studied in magnetized electron systems. [8j The lattice structure and its melting transition in 
charged-particle systems have been intensively studied. [9|, KL0| T he melting transitions from vortex lattices to vortex 



liquids have also been studied in the type II superconductors, jl lMl3| ] Vortex- lattice melting by quantum fluctuations in 
Bose-Einstein condensates on a one-dimensional optical lattice was theoretically studied. [14[ The Lindemann criterion 
of lattice vibration is often used whether the lattice is melted or not in many studies. 

On the other hand, point vortex systems have been intensively studied in fluid mechanics, especially as an ideal 
system of two-dimensional turbulence. The integrable, chaotic, and turbulent motions were found in some point vortex 
systems. [HI The statistical mechanics of point vortices have been studied since the pioneering work by Onsager. [l6l - [l8| 

In Bose-Einstein condensates, the vortex system can be approximated by a point vortex system if the repulsive 
interaction is sufficiently strong and the effect of the vortex core is negligible. In this paper, we discuss a melting 
transition of vortex lattices in general point vortex systems from the viewpoints of nonlinear dynamics and statistical 
mechanics. 

We perform some direct numerical simulations of point vortex systems and find that there is a transition from a 
layered distribution to a uniform distribution for the positional distribution of point vortices. We interpret the transi- 
tion as the melting transition and study the transition using a mean-field approximation for the layered distribution. 
A similar layered structure and the melting transition were found in a two-dimensional system of charged particles 
using Monte-Carlo simulations; [9] however, a melting transition in point vortex systems has not been numerically 
studied. The melting transition in a system of charged particles is somewhat similar to that in our model, because 
layered structures appear and the layered structures are destroyed at melting transitions. However, the dynamics of 
our point- vortex system is different from that in charged particles, in that the conservation law of angular momentum 
plays an important role in the point- vortex system. The equilibrium distribution is different from a simple canonical 
distribution owing to the additional conservation law. For example, an average rotational flow can exist even in an 
equilibrium state in the point vortex system, although such average flow is impossible in typical equilibrium states. We 
consider that the point vortex system is an interesting model of statistical mechanics in that long-range interactions 
and the additional conservation law are essential. 
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II. MODELS OF POINT VORTEX SYSTEMS 



We study point vortex systems because the equation of motion is simple, and a large-scale numerical simulation is 
possible in contrast to the numerical simulation of partial differential equations such as the Navier-Stokes and Gross- 
Pit aevskii equations. The point vortex systems are suitable for quantized vortices in superfluids because the circulation 
is quantized to a definite constant value, while the vorticity field changes continuously in most classical fluids. In 
three dimensions, a point vortex becomes a vortex filament. Vortex filaments were used to study vortex motions 
in superfluids and superfluid turbulences theoretically. [19| It is known that the approximation by point vortices and 
vortex filaments of vortices in superfluids is f airly good. Point vortices in two dimensions obey the equation of motion 
as explained in textbooks of fluid mechanics [20|: 

dxj _ 1 Kj(yj - yj) 



dt 2tt (pa - Xj) 2 + (yi - yj) 2 ' 

dyi ^ 1 v Kj{xi — Xj) 

~dt ~ ^2£( Xi - Xj )2 + ( yi - yj )2> [ > 

where i = 1, 2, • • • , TV, N is the total number of vortices, (xi,yi) denotes the two-dimensional position of the zth 
vortex, and K{ is the circulation of the zth vortex. Point vortices interact with other vortices through a long-range 
force equivalent to the two-dimensional Coulomb force; however, the equation of motion is not Newton's equation of 
motion as in a system of charged particles. Each point vortex just flows passively under the velocity field generated 
by other point vortices. For quantized vortices in Bose-Einstein condensates, the circulation K{ takes discrete values 
Ki = 27rnh/m, where n is an integer and m denotes the atomic mass. In most cases, n takes +1 or —1, because a 
vortex of large n tends to break up into n vortices of circulation 27vh/m. If the circulation k^s are all positive, an 
anticlockwise rotational flow of average angular velocity Q is induced on the average. The average rotational flow 
seems to disappear, if the system is observed from a rotational frame of angular velocity of — Q. In the rotational 
frame, the equation of motion is expressed as 

dxi _ J_ ^ K,(yi - yj) 



dt 2tt^ ( Xi - Xj) 2 + (yi - yj) 2 ' 

dt ~ 2K^.(x l -x 3 ) 2 + (y l -y 3 ) 2 ' [Z) 

Here, we have assumed that all the point vortices have the same circulation m = k = 2irh/m. If the spatial scale 
transformation x' = x/R^y' = y/R, where R is the length scale such as the radius of a container used to confine the 
vortices, and the temporal scale transformation t' = nt/(2TiR 2 ) are assumed, the equation of motion of point vortices 
is written in a dimensionless form as 

dxi yi- yj 



dt 1 f^( x i- x j) 2 + (Vi ~ Vj) 2 ' 



dyi 



^_^2 + ^_^2 > (3) 



dt 1 j^.(xi-Xj) 2 + (yi-yj) 2 ' 

where x' i: y' i: and t' are rewritten respectively as Xi,yi and t, and uo = 2i\R 2 ^tjn. The energy (Hamiltonian) E in this 
system is expressed as 

£ = 4£E ln ^ + E^ 2 > ( 4 ) 

i j^i i 



where = \fx\ + y 2 and Tij = y (xi — Xj) 2 + (yi — yj) 2 . (The energy is expressed as h 2 / (mR 2 )E if dimensions are 
recovered.) Equation (3) can be expressed using E as 

dxi _ dE dyi __ dE 
dt dyi ' dt dxi ' 

In this time evolution, the total energy E is conserved. The angular momentum L = ^2i{xi(dyi/ 'dt) — yi(dxi/dt)} 
is also conserved in the time evolution of eq. (3). The center of mass X = (1/N) ^2iXi and Y = (1/N) ^2iyi obeys 
dX/dt = uoY and dY/dt = -uoX. 
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In the approximation of point vortices, the effect of the vortex core is neglected. If the core size is sufficiently 
smaller than the average interval between neighboring vortices, the approximation might be good. The effect of the 
confinement by a harmonic potential as often used in Bose-Einstein condensates is not easily incorporated in point 
vortex systems. However, the confinement by a hard wall potential at r = R can be incorporated in point vortex 
systems by using antivortices at the mirror-image positions R 2 /r^, but for simplicity, we do not consider the effect in 
this paper. 

The point vortices including N vortices in a rotating container take a triangular lattice as a stable configuration. 
The energy E is minimized at the stable configuration. The stable configuration can be numerically obtained using 
the equations: 



dt ^ (xi - Xj) 2 + (yi - yj) 2 ' 
dyi = ( \ - Vi - Vj , s 

dt ^^2^ ( x ._ Xj) 2 + {y ._ yj) 2- W 

This is because eq. (6) is expressed as dxi/dt = —dE/dxi and dyi/dt = —dE/dyi, and the energy E decreases with 
time, and a state of the lowest energy is obtained as a stationary state in the time evolution of eq. (6). The position 
of the zth vortex at the lowest energy is expressed as (x^^yio)- Vortices are attracted to the origin (0,0) by the first 
term in eq. (6) and they interact with each other via the repulsive force expressed by the second term in eq. (6). As a 
result, the vortices find a stationary configuration such as a triangular lattice. Campbell and Ziff found various stable 
configurations of N point vortices. [2 ll l22j. Stationary solutions to eq. (6) are stationary solutions to eq. (3). Finding 
stationary solutions to eq. (3) through eq (6) is a useful method. 

The angular momentum L = J2 i {xi(dyi/dt) — yi(dxi/dt)} can be calculated as 

l = - $>(*? + y 2 ) } + £ ^-^fy^-y^ = _ uI + W^) 

i 7^ (xi - xj) 2 + ( Vi - Vj ) 2 2 

where / = ^ r 2 . At the stationary states, the quantity I satisfies / = N(N — 1)/(2cj), because L — 0. Assume that 
the point vortices take a triagular lattice at the stable configuration. If the interval between the nearest-neighbor 
vortices is denoted as the area of the elemental triangle of the lattice is expressed as S = V3d 2 /4. If the average 
radius of the entire triangular lattice is denoted as Ro and the average density of vortices is denoted as p, p is expressed 
as (3/6) /S = 2v // 3/(3<i 2 ), the total vortex number is N — ttRqP, and / is expressed as I = N(N — l)/(2u) — ttRqp/2. 
From these relations, typical length scales are evaluated as Ro = \/ (N — V)Juj and d = (27r/cj) 1//2 (l/3) 1 / 4 . 

If the point vortices are mixed up well by the high-dimensional chaos and there is no spatial correlation, we can 
apply the Debye-Hiickel theory of electrolyte solution [22| for the distribution function P(r), although there might 
be a criticism that a canonical distribution could not be applied to a small system of vortex number N = 7. In the 
Debye-Hiickel theory, the density P(r) depends only on the radius r from the center. The density of the electrolyte 
solution and the electric potential are determined self-consistently. That is, the electric potential is determined by 
the density of the electrolyte solution through the Poisson equation, and the density is determined by the canonical 
distribution under the electric potential. The density P(r) of the electrolyte solution around an electrode obeys the 
relation P(r) oc exp{— /30(r)}, where r is the distance from the electrode, (j) is the electric potential, j3 = q/(kBT), q 
is the charge of ions, ks is the Boltzmann constant, and T is the temperature. Similarly, the stationary distribution 
P(r) of point vortices is expected to obey a canonical distribution, 

P(r) = exp{-/30(r) - 7 (l/2)^r 2 }/Z, (7) 

where <j)(r) is an effective potential for the point vortex, and Z is determined from the normalization condition 
J °° P(r)2irrdr = N. The two parameters f3 and 7 are determined from the conservation laws of the interaction 
energy E\ = — | mr i,j an d the rotational energy E2 = (1/2)cj^ r 2 . If f3 = 7, P(r) has the form of the standard 
canonical distribution P(r) = ex.p(—/3E)/Z where E = E\ + E2. However, 7 takes a different value from j3 in our 
system. The distribution is a generalization of the distribution in the Debye-Hiickel theory in that the conservation 
law of angular momentum is taken into consideration. In our system, the interaction energy is specified and the 
effective inverse temperature j3 is determined from the conservation law of E\. The two quantities E\ and E2 are 
expressed in the mean-field approximation as 

poo 

E x = - / P(r)Q(r) (In r)2irrdr, 
Jo 

rOO 

E 2 = (1/2)cj / r 2 P(r)2irrdr. (8) 
Jo 
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FIG. 1: (a) Vortex configuration at the lowest energy state at E = Eq = —370.33 in a point- vortex system of 37 vortices. 
Plots of (xi(t),yi(t)) at t = 0.2 x n (n = 1,2, • • • ,750) at (b) E = -369.522, and (c) E = -361.055. (d) Distribution P(r) at 
AE = 0.329, 0.809 and 9.27. 

Here, we have assumed that P(r) is equal to the density of the point vortex, and Q(r) is defined as Q(r) = 
J Q r P(r')2iir' dr' . That is, E\ expresses the integration of the interaction energy of the point vortex at r with all 
the point vortices within a circular region of radius r. E2 expresses the rotational energy using the vortex density 
P(r). The effective potential (j) satisfies 

This equation is equivalent to the Poisson equation in two dimensions. It is because the logarithmic interaction among 
point vortices is equivalent to that among two-dimensional point charges, and the electric potential obeys the Poisson 
equation under the distribution of point charges. The quantity E2 is fixed to be N(N — l)/4, because I = ^2 i rf 
is equal to I = N(N — l)/(2uS) owing to the conservation law of the total angular momentum. The parameter 7 
is determined as E2 is equal to N(N — l)/4. The above equations can be numerically solved. If j3 = 00 or the 
temperature is zero, P{r) takes the form of a step function: P{r) = ujN/{tt(N — 1)} for r < ro = y{N — and 
P(r) = for r > ro- The length ro is the same as Rq previously evaluated for the triangular lattice in the lowest 
energy. This type of mean- field approximation was studied previously and reported in refs.17 and 18. 

We will study the dynamical and statistical behaviors of point vortex systems using eq. (3). Although the energy 
(Hamiltonian) (4) seems to depend on r explicitly, the system is effectively uniform, because the term (l/2)ur 2 is 
cancelled by the repulsive interaction among vortices. Assuming the uniformity, we use initial conditions randomly 
perturbed from the regular triangular lattice. That is, vortices are perturbed from the stable configuration (#io,2/io) 
as (xi(Q),yi(Q)) = (xio + r^^o + r^), where r X i and r y i are random numbers randomly chosen from a uniform 
probability distribution between — ro and ro- The standard deviation Aro = {XX r xi + ^i)/^} 1 / 2 is evaluated as 
Aro = ^2/3^0 • The parameter ro denotes the magnitude of the random perturbation. We have further adjusted the 
random numbers as £\ r xi = 5^ r yi = and ^{(xio + r xi ) 2 + (^0 + r yi) 2 } = N (N - l)/(2u). Under these initial 
conditions, the center of mass is always (0, 0) and the angular momentum is kept to be L = 0. The total energy E 
increases with ro- We have used these random initial conditions to change the total energy. As E is increased, chaotic 
dynamics is observed. The ergodicity and a thermal equilibrium state are expected, as in the molecular dynamics 
(MD) simulation using Newton's equation of motion. Our numerical simulation corresponds to the MD simulation 
of the lattice vibration of a two-dimensional triangular lattice where initial conditions are set to be in a randomly 
perturbed lattice structure. The temperature of the point vortex system increases with the internal energy E; however, 
no explicit form of the temperature is known in the point vortex system, while the temperature is expressed using 
the temporal average of kinetic energy in the MD simulation. We have performed numerical simulations of systems 
including various total numbers of vortices. We show a typical example of N = 37 in this paper to show the melting 
transitions from the viewpoint of dynamical systems including those with a relatively small total number of vortices. 

III. THIRTY-SEVEN VORTEX SYSTEM 

In this section, we study a point vortex system with thirty-seven vortices at uj = 3. Figure 1(a) displays a vortex 
configuration at the lowest energy state at E = £?o = —370.33. Figures 1(b) and 1(c) show plots of (xi(t),yi(t)) at 
t = 0.2 x n (n = 1, 2, • • • , 750) at (b) E = -369.522, and (c) E = -361.055. At E = -369.522, a four-layered structure 
appears. At E = —361.055, the layered structure is broken and a uniform distribution appears. Vortices come and 
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FIG. 2: Stationary solutions of the mean-field equation (10) at (a) f3 = 42.5, (b) = 7.5, (c) /3 = 5, and (d) /3 = 2.5. 

go between the outer three layers for AE = E — Eq > 0.301. Figure 1(d) shows the distribution P(r) of the vortices 
in the radial direction at AE = 0.329, 0.809 and 9.27 by the numerical simulation. The peak structures of P(r) at 
AE = 0.329 and 0.809 correspond to the layered structures, and the flat distribution at AE = 9.27 corresponds to 
the liquidlike state. 

The melting transition can be treated qualitatively by the mean-field approximation even for this four-layered 
structure. The probability distributions for the four layers are expressed as Pi, P2, P2, and P4, and the potentials are 
denoted as 01,02,03, and ^4. 

P t (r) = exp{-^(r)- 7 (l/2)^r 2 }/Z z , 

E x = - / {P(r)Q(r)-Vg < P i (r)Q < (r)}(liir)27rrdr, 
J ° i=i 

PCX) 

E 2 = (l/2)o; / r 2 P(r)2irrdr, 
Jo 

V ^ = -rl{ rd t)=-^ P -^ (10) 

where Z^s are normalization constants for P^, Qi = J Q r Pi(r / )2irr / dr\ P(r) = Xlt=i^( r )' anc ^ Q( r ) = Y^t=iQ i ( r )- 
The quantities ^'s denote the proportion of self-interaction in each layer: q\ — 1,^2 = 1/6, <73 = 1/12, and = 1/18. 
Stationary solutions to eq. (10) were numerically obtained. Figure 2(a) shows the four-layered solution for f3 = 42.5 > 
f)\. The interval between the layers is about 0.95, which is comparable to the result of the direct numerical simulation. 
The layered structure disappears from the outer layers. The outer two layers are merged into one layer in Fig. 2(b) at 
j3 = 10.3 < P\. The outer three layers are merged into one layer in Fig. 2(c) at j3 = 5.9 < fa. The layered structure 
disappears completely and changes into a liquidlike state at f3 = 2.5 < ^3. 

IV. SUMMARY 

We have performed a direct numerical simulation of point vortex systems. We have found that layered structures 
appear when the energy increases. The layered structures are gradually broken and a uniform distribution for the 
position of vortices is realized when the energy increases sufficiently. Taking the layered structure into consideration, 
we have derived a mean- field equation for the distribution P^ for each layer. The distribution is assumed to obey 
a canonical distribution based on the two conservation laws of the total energy and the total angular momentum. 
The interval between the neighboring layers can be evaluated in the mean- field theory for the layered structure. The 
melting transition from a layered structure to a liquid state appears naturally in the mean-field approximation. 
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